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i^^' Abstract. Various classes of stable finite difference schemes can be con- 

' structed to obtain a numerical solution. It is important to select among 

all stable schemes such a scheme that is optimal in terms of certain ad- 
ditional criteria. In this study, we use a simple boundary value problem 
for a one-dimensional parabolic equation to discuss the selection of an 
approximation with respect to time. We consider the pure diffusion equa- 
tion, the pure convective transport equation and combined convection- 
■^C . diffusion phenomena. Requirements for the unconditionally stable finite 

h^ ' difference schemes are formulated that are related to retaining the main 

^^ I features of the differential problem. The concept of SM stable finite dif- 

C/3 . ference scheme is introduced. The starting point are difference schemes 

constructed on the basis of the various Fade approximations. 
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When time-dependent problems of mathematical physics are solved numerically, 
\f^ • much emphasis is placed on computational algorithms of higher orders of ac- 

curacy (e.g., see [ll2j ). Along with improving the approximation accuracy with 
p^ • respect to space, improving the approximation accuracy with respect to time is 

(^ [ also of interest. In this respect, the results concerning the numerical methods 

for ordinary differential equations (ODEs) |3l4j provide an example. Taking into 
account the specific features of time-dependent problems for PDEs, we are inter- 
. ■ ested in numerical methods for solving the Cauchy problem in the case of stiff 

rS I equations |5I6I7J . 

jrt ' When time-dependent problems are solved approximately, the accuracy can 

be improved in various ways. In the case of two-level schemes (the solution at 
two adjacent time levels is involved), polynomial approximations of the scheme 
operators on the solutions are used explicitly or implicitly. The most popular rep- 
resentatives of such schemes are Runge-Kutta methods [718] , which are widely 
used in modern computations. The main feature of the multilevel schemes (mul- 
tistep methods) manifests itself in the approximation of time derivatives with a 
higher accuracy on a multipoint stencil. A characteristic example is provided by 
multistcp methods based on backward numerical differentiation [9]. 

Various classes of stable finite difference schemes can be constructed to obtain 
a numerical solution [lOlllj . It is important to select among all stable schemes 
such a scheme that is optimal in terms of certain additional criteria. In the theory 



of finite difference schemes, there is the class of asymptotically stable schemes 
(see |f2lf3] ) that ensure the correct long-time behavior of the approximate so- 
lution. In the theory of numerical methods for ODEs (see J7I9| ). the concept 
of L-stability is used, which reflects the long-time asymptotic behavior of the 
approximate solution from a different point of view. 

In [13] there arc considered properties of two-level difference schemes of high 
order approximation for the approximate solution of the Cauchy problem for 
evolutionary equations with self-adjoint operators. The simplest boundary value 
problem for the one-dimensional parabolic equation serves as a basic problem. 
The concept of SM stability (Spectral Mimetic stability) of a difference scheme is 
introduced. This property is connected with the behavior of individual harmonics 
of the approximate solutions. 

In this paper, we continue to study the SM properties of difference schemes for 
the approximate solution of unsteady problems of mathematical physics. For the 
model boundary value problem for one-dimensional parabolic equation there are 
considered the spectral characteristics of approximations in space and in time. 
In particular, good approximation properties are observed for the convection 
operator of third order in space. Features two-level schemes of higher order of 
approximation in time, based on the Fade approximation, are considered for 
solving problems of mathematical physics with symmetric and skew-symmetric 
operators. 

2 Problem formulation 

We consider finite-dimensional real Hilbert space H, where the scalar product 
and the norm is (•, •) and || • ||, respectively. Let u{t) (0 < t < T > 0) is defined 
as the solution of the Cauchy problem for evolutionary equation of first order: 

^+Au = f{t), 0<t<T, (1) 

m(0) = uo. (2) 

The right part f{t) G H oi equation ([T]) is given and yl is a linear non-negative, 
in general, nonself-adjoint operator from H in H {A = A{t) > 0) depending on 
t. 

For problem (IT|), ^ the estimate of stability is easily established. Taking 
into account the skew-symmetric property of operator A, we have the equality 

II ii^ll"ll ff \ 

M^^ = iLu). 

With regard 

(/,^)<lh^ll 11/11 

we obtain a simple estimate of stability for the solution of ([T|) , ([2]) with respect 
to the initial data and right hand side: 

\Ht)\\ < hoW + I \\fm\de. (3) 



We would like to preserve these properties of the differential problem after tran- 
sition to a discrete analogue of problem ([1]) , ([2]) . 

The main attention in our discussion is given to unsteady boundary value 
problem for partial differential equations. In this context, we can associate the 
Cauchy problem ([1]), ([2]) with application of the method of lines (approximation 
in space). Given the importance for applications, wc will conduct our consid- 
eration on an example of a boundary value problem for the one-dimensional 
parabolic equation of second order. Let a sufficiently smooth function u{x,t) 
satisfies the equation 

-:^ + Cu = 0, < x < 1, < i < T (4) 

ot 

and the initial conditions 

u{x,0) = uo{x), < a; < 1. (5) 

The periodicity of the spatial variable is assumed: 

u{x + l,t)=u{x,t), 0<t<T. (6) 

We associate operator C with the convection-diffusion equation, defining 

£ = XC + (1 - x)P (7) 

for some constant < x ^ 1- Here the operators of convective and diffusive 
transport are defined according 



C"=-' 2?^. = --J. (8) 



du d'^u 

dx ' dx^ 

If X = 1 equation (J4|) is the convection transport equation whereas if x = it is 
the diffusion equation. 

The discrete problem should inherit the main properties of the differential 
problem. In model problem (|4])-(|6]) the skew-symmetric property of operator C 
as well as the self-adjoint and non-negative properties of operator D should be 
preserved 

C = -C*, V = V*>0 (9) 

in space £2(0, 1) for functions satisfying ©. Stability of the solution of the corre- 
sponding problem ([T|). ^ (estimate (|31)) is provided by similar properties of the 
grid analogs of convective and diffusive transport operators. In our research, we 
are concentrating on the spectral characteristic of the solution (Spectral Mimetic 
Properties for grid approximations) , when considering the behavior of individual 
harmonics of the approximate solution. 

3 SM properties of the approximation in space 

Let us introduce a uniform grid with step h: 

Xi = ih, i = 0, ±1, ±2, ..., Mh = 1, 



uj^{x, I 1 = 0,1,..., M-1}. 

We use standard index- free notations of the theory of difference schemes [TO] . 
Let w = Wi = w{xi), and for the left, right and central difference derivatives we 
set 

= j^ , o+w = , 

9ow = 2(<9-^ + 9+w) = ^^^^ 

respectively. 

After approximation in space we put for problem (|4])-(|6l) the corresponding 
discrete problem 

-Ji+Ay = 0, xeu, o<t<r, (10) 

at 

y{x,0) ^uo{x), xeuj. (11) 

Some key possibilities in the choice of grid operator A with the focus on the 
properties of the differential operator C should be noted. 

We define Hilbert space H = L2{oj) of periodic grid functions {y{x + 1) = 
y{x)) with the inner product and the norm 



{y,w)^^y{x)w{x)h, ||j/f = (y,y). 



To guarantee stability of solution (fTO|). ([TT]) . operator A must be non- negative 
{A > 0) in H. Conservatism (neutral stability) communicates directly with the 
skew-symmetric property of operator yl (yl = — ^4*) in H. 

For the convection equation (x = 1, -C = C) with / = the norm of the 
solution of problem (g])-® does not change in time: 

||u(i)|| = ||,.o||. (12) 

Equality P^ reflects the conservation property of the solution (conservation 
law), the neutral stability of the solution. 

1. Upwind (directional) approximations of first order. To approximate the 
convective terms (see, e.g., [1115116] ). the upwind first-order approximations are 
traditionally widely used. In this case, the grid convection operator C has the 
form 

C = d-. (13) 

Operator C defined in ([T3|) is non- negative (C > 0). In this case, for the solution 
of problem (fTO|) . (fTT|) the estimate 

||j/||<ho||, 0<t<T (14) 

is true. 

2. Central- difference approximation. Another well-known variant is to use 
approximations of second order where 

C^do. (15) 



In this case we have C = -C* and for problem (dU]), (HH, dH]) holds 

\\y\\ = \\uol 0<t<T. (16) 

3. Upwind second- order approximations. When choosing approximations of 
higher order (second and above) for the convective terms, we are trying at least 
partially to preserve the properties of the first order approximations, which are 
connected primarily with the monotonicity (fulfillment of the maximum prin- 
ciple) . The most interesting attempts in the class of linear approximations are 
associated with the use of approximations with the upwind differences of second 
order |2I17| . For our problem (|3])-(jn]), we have 

^^^ 2h ■ 

Using previously introduced operator notations we obtain 

C = 9_ + |a_a_. (17) 

Operator C > and so for the solution of problem ([TU)) , ([TT|) estimate P^ holds 
again. 

4- Approximations of third order. In computing practice third order approxi- 
mations are not in common use. In fact, they are only mentioned (see, e.g., j2|4j ) 
without any meaningful analysis. In this case the difference convection operator 
can be written in the form 

C = do-^d.d.d+. (18) 



In index notation equation pSI) takes the form 

'iVi+i + 3yi - 6yi-i + yi-2 
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Operator C > and its energy (equal to {Cy,y)) is three times less energy 
operator C, which is defined by rule (|17p (upwind second-order approximations). 

The stability conditions (neutral stability) of the considered approximations 
of convection transfer are associated with the general properties of the operator 
(non- negativity, skew-symmetric property). More detailed information gives us 
the spectrum of the difference operator, its proximity to the spectrum of differ- 
ential operator. This inheriting the properties of the differential problem for the 
difference problem at the spectral level we associate [TJ| with SM properties. 

Consider the corresponding differential problem for eigenvalues and eigen- 
functions. For operator C we have 

dv 

— = Xv, 0<a;<l, (19) 

dx 

v{x + l)=v{x). (20) 



The solution of spectral problem p^ .ipH )) is 

Am = i2Trm, 

v^{x) = e''^^\ m = 0,±l,±2,.... 
For solution of problem ([I])-® we obtain the representation 

u{x,t)= ^ (uo,Wm)e~^™*Wm(a;), (21) 

m— — oo 

where 

{uo,V,n)= Uo{x)v,n{x)dx, 771 = 0, ±1, ±2, ... 

Jo 
are the coefficients of expansion for function uq (x) . 

We now consider the appropriate spectral problems of the grid problem 

Cv = nv (22) 

with the above-mentioned approximations of convcctive transport. For dcfinite- 
ness, we assume that M is odd. Eigcnfunctions of problem (|22p for the considered 
difference operators ([T^]) . ((TS|) . p7|) and ((T5)) have the form 

ii;m(s) = e'2™-, xecu, m = 0,±l,±2,...,^^—^. (23) 



For difference operator C, defined according to p3)) . the eigenvalues have the 
form 

fJ-m^ j^ , m = 0,±l,±2,..., — - — . 

For the imaginary and real parts we obtain 

R-C fim — — sin (TTTTiAi), Im fi,fn = — sin(27r?Ti/i), m = 0, ±1, ±2, ..., . (24) 

h h 2 



For central-difference approximations (JTSj) we obtain 
Re fi„i ~ 0, Im /i„j = — sin(27r?7i/i), to = 0, ±1, ±2, ..., . (25) 

Comparing p4p and (|25p shows that the imaginary components of the spectrum 
of central-difference approximations and upwind difference approximations co- 
incide. For upwind approximations we have positive real parts of the spectrum, 
which cause the dissipative properties of such approximations. 

Dissipative properties demonstrate also approximations pT|) . (|18p . For the 
upwind second-order approximations we have 

Re fim = —icos(2-Kmh) — 1)^, Im /!„ = — sa\{2-Kmh)(2 — cos(27r?Ti/i)) (26) 
h h 
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with the above-mentioned values of m. For the third order approximations it is 
easy to obtain 



Re /i,„ = — (cos(27rTO/i) 



1 



Im ^r. 



1 



sin(27rm/i)(4-cos(27rm/i)) (27) 



respectively. 

An illustration of the spectrum of grid convection operator ([M)) -(l?7 |) is shown 
in Fig lll2l for M = 31. In particular, the main disadvantage of the scheme with 
directional differences of first order is associated with substantial dissipation of 
low harmonics whereas dissipative properties of scheme (|17p , (|18p are connected 
primarily with high harmonics. Approximation of third-order for (|18p is relatively 
well reflects the spectral properties of the differential problem. Its dissipative 
properties work only for high harmonics and is weak for the most important low 
harmonics of the difference solution. 

A similar analysis has been performed for the diffusion equation, where x = 1, 
£ = 2?. In this case, we investigated the non-negativity and self-adjointness of 
discrete diffusion operator C and its spectral properties. 

1. Approximation of second order. The standard approximation at the three- 
point stencil leads us to 

D = -d+d-. (28) 



2. Approximation of fourth order. At the extended stencil we can use 



D 



-d+d- 



h^ 
12' 



d+d-d+d-. 



(29) 
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Using dlHl), dMl), we have D ^ D* > in H. 

The spectrum of these operators D is real with the same eigenfunctions that 
for C . For the eigenvalues we have 



/?^^" -M 



TO ==0,1,..., A/- 1 



if (EHl) and 



^'- = V^'''' -M 



l + -sin' 



mix 



TO 



0, 1,...,M-1 



if approximation ([^^ are selected. For the differential operator we have A^ = 
47r^TO^, m ~ 0,1,.... As expected, approximation (P^ gives us the better ap- 
proximations for the spectrum of the differential operator of diffusive transport. 



4 SM properties of the approximation in time 

We'll use the two-level difference schemes for the approximate solution of (fTO]) . 
(jlip . Define a uniform grid in time with the time-step r 

aJ^ = cj^ U {T} = {t„ = nr, n = 0, 1, ..., A^, tN^T] 

and let ijn = y{tn), tn — nr. For the exact solution of problem pH]) . PT|) . at 
transition from time level t„ to new time level t„+i we have 

(M-l)/2 

y{x,tn+i) ^ e~^''y{x,tn) ^ ^ (y(a;, i„), w„)e"^"^w„(a;). (30) 

m=-[M-l)/2 



Two- level difference scheme for problem pH)) , pT|) is written in the canonical 
operator-difference form 

g ^"+'~^" +Ay„ = 0, n = 0,l,... (31) 

r 

with some operators A and B. In the Samarskii theory of stability of operator- 
difference schemes [lOlllll^ stability conditions in the various norms are for- 
mulated in the form of operator inequalities for the A, B. 
Difference scheme ([51]) is written as follows 

Vn+i = Syn, n = 0,l,..., (32) 

where 

S ^E-tB-^A (33) 

is the operator of transition from one time level to another level, which, in 
general, may depend on n. 

We restrict ourselves to the simplest difference approximation in time for 
problem ([T0|. pTj) . which lead to the transition operator 

S^s{tA), (34) 

where s{z) is a function of stability |&17J. With constraints (|M)) {tA = {tA){tA), 
B = B(tA)) the stability conditions in Hilbcrt spaces are easily verified on the 
basis of only properties of function s{z). Let A > SE, then 

\\s{tA)\\ < max |s(z)|. 

Re z>5r 

and self-adjointness of operator A is not assumed. 

In the case of (p4)) for the approximate solution at the new time level we have 
the representation 

{M-l)/2 
y{x,tn+l)= ^ {y{x,tn),Wm)s{,lJLmT)Wm{x). (35) 

m=-(A/-l)/2 



Quality of difference approximations in time is estimated by comparing (j35p with 
representation ([5(I| for model problem ([TU| . ([TT|) . The comparison is performed 
at the level of behavior of individual harmonics and so we are talking about the 
SM properties for approximation in time. 

For convection problem (x = 1, i2 = C) the spectrum is purely imaginary, and 
the solution is neutrally stable. After approximation in space using the above 
directional differences a typical situation is where the imaginary part of the 
spectrum is complemented by real part. When choosing approximations in time 
for the considered problems with skew-symmetric operators, wc must monitor 
the behavior of the main imaginary part of the spectrum. This means that in 
problem ^, ^ 

A = A^ + A^, Ao = Al = ]^{A + A*), A, ^ -Al = ]^{A- A*) 



operator Ai is main in the sense that Aqi/ — > as /i — ;• for sufficiently smooth 
y. The supporting real part of the spectrum associated with operator Aq, is gen- 
erated by the approximations in space and plays a minor role in these problems 
(it should not lead to instability of the difference solution) . 

The difference scheme for convection problem pO)) . (fTTj) . in which operator 
yl > and its antisymmetric part has the major role, is called as the SM stable 
if the difference scheme is stable and neutrally stable a.t A — ~A*. 

Two-level difference schemes of higher order accuracy for time-dependent 
linear problems we will construct on the basis of the Pade approximations for 
the operator (matrix) exponent e~ . For e~^ we have 

e-^ = i?;,„(z) + 0(z'+™+i), RUz) ^ ^4^, 

Qlm\Z) 

where Pim{z) and Qimiz) arc polynomials of degree I and m, respectively: 

l\ ' {l + m-k)\ , 

^'"^^~ (/ + m)!^ k\{m-k)l ■ 

For equation (|fOp the application of Fade approximations corresponds to the 
two-level scheme 

QiUrA f"^^ ~ ^" + -{QhnirA) - P;,„(t^));/„ = 0, n = Q,l,.... (36) 

T T 

In canonical form (|3fp difference scheme ([55]) corresponds to the choice 

A=-{Qirn{TA)-Pi„,iTA)), B = Qu„{tA). (37) 

r 

Difference schemes for problem ([TU)) , pT|) with yl > on the basis of Fade 
approximations are stable (absolutely stable) (estimate ||yn-i-il| < ||yn|| holds) at 
I < m [617] . It is only necessary to highlight among such schemes the SM-stable 
schemes. 

In the simplest case m = 1 we have 

Roiiz) = -^ = e-^ + Oiz'), 
1 + z 

Rii{z)^\-i^=e-^ + 0{z^). 

1 + ^z 

Approximations Roi{z) corresponds to application for the approximate solution 
of problem (fTO|) , (fTTj) the purely implicit scheme 

^"+^"^" +y4y„+i^0, n = 0,l,.... (38) 



The application of the symmetric scheme (Crank-Nicholson) 

?/»+l-y» ^^ j/n+l+yn ^Q^ „ = 0,1,.... (39) 

r 2 

corresponds to the choice of approximation Rii{z). 

The condition of neutral stability of ||?/„+i|| = \\yn\\ for two-level scheme p2p 
will be satisfied at \\S\\ = 1. Taking into account ([M| and for A = —A* this 
corresponds to the case 

Re z = 0. (40) 
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we have 




|i?0] 


-Wl = 


1 




v/i + y' 



Thus, the condition of neutral stability is not satisfied — the purely implicit is 
not SM stable for problems with the main skew-symmetric operator. While for 
symmetric scheme pop we obtain 

|i?ii(z)| = l, z = iy. 

Thus, this scheme is SM stable for the investigated class of problems. 

We can make similar conclusions for schemes with Fade approximations at 
m > 1. Only a scheme that is based on approximation Rmm is SM stable for 
problems with the main skew-symmetric operator. At using Fade approxima- 
tions with / < m the scheme demonstrates dissipative properties due to the 
approximation in time. Only at I = m the corresponding scheme is neutrally 
stable. 

A similar analysis is carried out fsee|14j) for the diffusion equation. In prob- 
lem pUj) . ([TT]) with yl = yl* > amplitudes of harmonics with higher numbers 
damp more quickly in compare with amplitudes of harmonics with lower num- 
bers (spectral monotonicity) and damp to zero as t — ^ oo (asymptotic stability). 
Such a behavior of the approximate solutions we associate with the SM prop- 
erties of approximation in time for the solution of problems with self-adjoint 
operators. 

We assume that the difference scheme for problem ((T0| , (fTTj) with yl = yl* > 
is SM stable if it is spectrally monotonic and asymptotically stable. 

Difference schemes based on the Fade approximation i?i,„ are SM stable at 
I = 0. Furely implicit scheme ([55)) belongs to this class of schemes, whereas the 
symmetric scheme is conditionally SM stable. 

The main conclusion of the study is that the approximate solution of prob- 
lems with skew-symmetric operators we must use such approximations in time, 
which are based on the Fade approximations Rmm{z)- For problems with self- 
adjoint operators it is more preferred to use Fade approximation i?om(^)- For 
problems with general non-selfadjoint operators, approximation in time can be 
constructed via decomposition into the self-adjoint and skew- symmetric com- 
ponents and further constructing different approximations for them, based on 
special splitting schemes. 
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